!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine heat_conduction(ro, mx, my, mz, bx, by, bz, en, dx, dy, dt, nx, ny)
!==================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx, ny

   real(8), intent(in) :: dx, dy, dt
   real(8), dimension(nx, ny), intent(inout) :: ro, mx, my, mz, bx, by, bz, en

   integer(4) :: i, j !, ip1, jp1

   real(8) :: s_en, dx2, dy2
   real(8), dimension(nx, ny) :: en_pre, hc_x, hc_y

!----------------------------------------------------------------------------------------------------------------------------------|
   dx2 = 2.0d0 * dx
   dy2 = 2.0d0 * dy

!----------------------------------------------------------------------------------------------------------------------------------|
!  first step
!----------------------------------------------------------------------------------------------------------------------------------|
   call model_hc(hc_x, hc_y, ro, mx, my, mz, bx, by, bz, en, dx, dy, nx, ny)

!   do j = 1, ny
!   do i = 1, nx
!      tmp = (mx(i, j) * mx(i, j) + my(i, j) * my(i, j) + mz(i, j) * mz(i, j)) / ro(i, j) / 2.0d0 +                                 &
!            (bx(i, j) * bx(i, j) + by(i, j) * by(i, j) + bz(i, j) * bz(i, j)) / 2.0d0
!      te(i, j) = (en(i, j) - tmp) * gmm1 / ro(i, j)
!      if (te(i, j) .le. 0.0d0) then
!         te(i, j) = eps
!         en(i, j) = te(i, j) * ro(i, j) / gmm1 + tmp
!      endif
!      hc_x(i, j) = 0.0d0
!      hc_y(i, j) = 0.0d0
!   enddo
!   enddo

!   do j = 2, ny - 1
!   do i = 2, nx - 1
!      b2 = bx(i, j) * bx(i, j) + by(i, j) * by(i, j)
!      if (b2 .lt. eps) then
!         hc_x(i, j) = min((te(i + 1, j) - te(i - 1, j)) / dx2, 1000d0) * hc0 * te(i, j) ** (5.0d0 / 2.0d0)
!         hc_y(i, j) = min((te(i, j + 1) - te(i, j - 1)) / dy2, 1000d0) * hc0 * te(i, j) ** (5.0d0 / 2.0d0)
!      else
!         te_grad_b = min((te(i + 1, j) - te(i - 1, j)) / dx2, 1000d0) * bx(i, j) +                                                 &
!            min((te(i, j + 1) - te(i, j - 1)) / dy2, 1000d0) * by(i, j)
!         hc_x(i, j) = te_grad_b * bx(i, j) / b2 * hc0 * te(i, j) ** (5.0d0 / 2.0d0)
!         hc_y(i, j) = te_grad_b * by(i, j) / b2 * hc0 * te(i, j) ** (5.0d0 / 2.0d0)
!      endif
!   enddo
!   enddo

   do j = 2, ny - 1
   do i = 2, nx - 1
         s_en = (hc_x(i + 1, j) - hc_x(i - 1, j)) / dx2 + (hc_y(i, j + 1) - hc_y(i, j - 1)) / dy2
         en_pre(i, j) = en(i, j) + s_en * dt
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  second step
!----------------------------------------------------------------------------------------------------------------------------------|
   call model_hc(hc_x, hc_y, ro, mx, my, mz, bx, by, bz, en_pre, dx, dy, nx, ny)

!   do j = 1, ny
!   do i = 1, nx
!      tmp = (mx(i, j) * mx(i, j) + my(i, j) * my(i, j) + mz(i, j) * mz(i, j)) / ro(i, j) / 2.0d0 +                                 &
!            (bx(i, j) * bx(i, j) + by(i, j) * by(i, j) + bz(i, j) * bz(i, j)) / 2.0d0
!      te(i, j) = (en_pre(i, j) - tmp) * gmm1 / ro(i, j)
!      if (te(i, j) .le. 0.0d0) then
!         te(i, j) = eps
!         en_pre(i, j) = te(i, j) * ro(i, j) / gmm1 + tmp
!      endif
!      hc_x(i, j) = 0.0d0
!      hc_y(i, j) = 0.0d0
!   enddo
!   enddo

!   do j = 2, ny - 1
!   do i = 2, nx - 1
!      b2 = bx(i, j) * bx(i, j) + by(i, j) * by(i, j)
!      if (b2 .lt. eps) then
!         hc_x(i, j) = min((te(i + 1, j) - te(i - 1, j)) / dx2, 1000d0) * hc0 * te(i, j) ** (5.0d0 / 2.0d0)
!         hc_y(i, j) = min((te(i, j + 1) - te(i, j - 1)) / dy2, 1000d0) * hc0 * te(i, j) ** (5.0d0 / 2.0d0)
!      else
!         te_grad_b = min((te(i + 1, j) - te(i - 1, j)) / dx2, 1000d0) * bx(i, j) +                                                 &
!            min((te(i, j + 1) - te(i, j - 1)) / dy2, 1000d0) * by(i, j)
!         hc_x(i, j) = te_grad_b * bx(i, j) / b2 * hc0 * te(i, j) ** (5.0d0 / 2.0d0)
!         hc_y(i, j) = te_grad_b * by(i, j) / b2 * hc0 * te(i, j) ** (5.0d0 / 2.0d0)
!      endif
!   enddo
!   enddo

   do j = 2, ny - 1
   do i = 2, nx - 1
         s_en = (hc_x(i + 1, j) - hc_x(i - 1, j)) / dx2 + (hc_y(i, j + 1) - hc_y(i, j - 1)) / dy2
         en(i, j) = (en(i, j) + en_pre(i, j) + s_en * dt) * 0.5d0
   enddo
   enddo

!   do j = 1, ny - 1
!      jp1 = j + 1
!      do i = 1, nx - 1
!         ip1 = i + 1
!         b2 = bx(i, j) * bx(i, j) + by(i, j) * by(i, j)
!         if (b2 .lt. eps) then
!            hc_x(i, j) = min((te(ip1, j) - te(i, j)) / dx, 1000d0) * hc0 *                                           &
!               ((te(ip1, j) + te(i, j)) / 2.0d0) ** (5.0d0 / 2.0d0)
!            hc_y(i, j) = min((te(i, jp1) - te(i, j)) / dy, 1000d0) * hc0 *                                           &
!               ((te(i, jp1) + te(i, j)) / 2.0d0) ** (5.0d0 / 2.0d0)
!         else
!            te_grad_b = min((te(i + 1, j) - te(i - 1, j)) / dx2, 1000d0) * bx(i, j) +                                             &
!               min((te(i, j + 1) - te(i, j - 1)) / dy2, 1000d0) * by(i, j)
!            hc_x(i, j) = te_grad_b * bx(i, j) / b2 * hc0 * te(i, j) ** (5.0d0 / 2.0d0)
!            hc_y(i, j) = te_grad_b * by(i, j) / b2 * hc0 * te(i, j) ** (5.0d0 / 2.0d0)
!         endif
!      enddo
!   enddo

!   do j = 2, ny
!   do i = 2, nx
!         s_en = (hc_x(i, j) - hc_x(i - 1, j)) / dx + (hc_y(i, j) - hc_y(i, j - 1)) / dy
!         en(i, j) = en(i, j) + s_en * dt
!   enddo
!   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine heat_conduction
